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Abstract 

We analyze the dynamics of a class of cosmological solutions of the Einstein- Vlasov equa- 
tions. These equations describe an ensemble of collisionless particles (which represent galaxies 
or clusters of galaxies) that interact gravitatively through Einstein's equations of general rel- 
ativity. The cosmological models we consider are spatially homogeneous, of Bianchi type IX, 
and locally rotationally symmetric (LRS). We prove that generic solutions exhibit an oscilla- 
tory approach toward the singularities (the 'big bang' in the past and the 'big crunch' in the 
future) ; this is in contrast to the behavior of Einstein- vacuum or Einstein-Euler solutions. To 
establish this result we make use of dynamical systems theory; we introduce dimensionless dy- 
namical variables that are defined on a compact state space; in this formulation the oscillatory 
behavior of generic solutions is represented by an approach to heteroclinic cycles. 



1 Introduction 



A large number of recent rigorous results on the dynamics of cosmological solutions (i.e., models 
for the 'universe') of the Einstein field equations is due to the successful application of the theory of 
dynamical systems to general relativity. The Einstein equations form a highly non-linear system of 
partial differential equations that describe the evolution of the metric of the space-time, the latter 
being represented by a four-dimensional Lorentzian manifold. In the case of spatially homogeneous 
solutions, which are of interest in cosmology, the Einstein equations reduce to a system of non-linear 
ordinary differential equations, which can be analyzed using the powerful methods of dynamical 
systems theory. This is typically achieved by going over from the standard metric variables of 
the Einstein equations to dimensionless variables via conformal rescalings and the introduction 
of normalizations to regularize the equations and to obtain an autonomous finite-dimensional 
dynamical system over a compact state space. We refer to [3T] for an overview on the theory of 
dynamical systems in cosmology. 
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This paper concerns cosmological solutions of the Einstein equations coupled to collisionless matter 
(Vlasov matter). Solutions of the Einstein- Vlasov equations represent ensembles of massive par- 
ticles (like stars in a galaxy) that interact through the gravitational field they create collectively. 
In cosmological applications, the particles are thought to represent galaxies (or galaxy clusters) in 
the universe. 

According to the standard theory for the evolution of the universe, massless particles (photons) 
account for most of the energy density in the universe during the time of 'radiation dominance', 
which begins at about 1 second after the big bang and ends at the time of decoupling between 
radiation and matter about 10 5 years after the big bang. From this time onward, the universe is 
'matter dominated' and galaxies begin to form about 10 6 years after the big bang. It is immediate 
that the collisionless matter model that describes ensembles of massive particles is a useful tool 
to model structure formation. The time of 'radiation dominance', on the other hand, suggests 
the study of ensembles of massless particles (where in the approach to the big bang singularity 
collisions are of course expected to play an increasing role — this is not captured by the collisionless 
matter model). Interestingly enough, the dynamics toward the big bang singularity in the case of 
massive particles is qualitatively the same as the dynamics of the massless particles model. This 
seems to be a generic property of Vlasov matter, see [TT1 [20] . 

The existence of a big bang (an 'initial singularity') is predicted by the singularity theorems 
of general relativity |10) . Under very general assumptions, there will exist causal curves that 
'terminate' at a singularity, and physical quantities representing curvature or the energy density 
of the matter will diverge along these inextendible curves. The detailed characterization of these 
singularities is an important problem in general relativity and cosmology. One interesting question 
in this respect concerns the details of the divergence of the curvature toward the singularity (e.g., 
along a distinguished congruence of curves representing the (spacetime) trajectories of the matter). 
The well-known BKL conjecture [3] states that the approach to the singularity will in general be 
'oscillatory', which means that appropriately rescaled (curvature) quantities will oscillate (with 
a rapidly increasing frequency) instead of converging monotonically toward the singularity. The 
paradigm of this type of behavior is the so-called 'Mixmaster' behavior that originates from the 
study of spatially homogeneous spacetimes. We refer to sections 5 and 6 of the textbook [5TJ for 
a good overview. Another important problem concerns the question of whether the asymptotic 
dynamics of solutions toward the singularity are sensitive to the choice of matter model or not 
(i.e., whether "matter matters"). Although for matter models like perfect fluids the latter seems 
to be the case, there are matter models such that the solutions of the Einstein-matter equations 
exhibit an asymptotic behavior that is different from that of vacuum solutions. The collisionless 
(Vlasov) matter model is a prime example, and the results of the present work are results in this 
vein. 

To obtain a detailed characterization of the dynamics of solutions of the Einstein- Vlasov equa- 
tions, it is necessary to assume a high degree of symmetry as, e.g., spatial homogeneity. The 
global dynamics of spatially homogeneous solutions of the Einstein- Vlasov system has been stud- 
ied extensively, see, e.g., [3 [TTJ [TBI HZl UHl [20] for applications of dynamical system theory to the 
problem. The reformulation of the Einstein (-Vlasov) equations as a dynamical system has proved 
to be highly advantageous. The reason for the success of dynamical systems methods lies in the 
fact that the behavior of the spacetime geometry in the neighborhood of an (initial) singularity is 
determined by the nature of the a-limit set of the dynamical system; see section 5.3 of PL for a 
good introduction. In the simplest cases, the a-limit set is an isolated fixed point; typically, this 
corresponds to the spacetime curvature growing monotonically. If the a-limit set is more com- 
plicated, however, for example a heteroclinic cycle with a finite number of fixed points, then the 
approach to the singularity will be 'oscillatory'. The paradigm of oscillatory behavior, the Mix- 
master behavior described in the BKL conjecture, is even more intricate — it is induced by infinite 
heteroclinic chains. 
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The present work extends the results of [19] on massless Vlasov matter and [20] on massive particles. 
The results of these references concern the dynamics of Einstein- Vlasov solutions that satisfy the 
most restrictive symmetry assumptions (LRS Bianchi type I, II and III); the analysis is based 
on techniques from dynamical systems theory in conjunction with the use of Hubblc-normalized 
dimensionless variables, see [21]. In the present paper we refine these techniques to investigate a 
class of solutions exhibiting a larger number of (true) degrees of freedom: LRS Bianchi type IX; we 
refer to section [5] for a definition. Our analysis employs a different set of dimensionless variables 
to regularize the equations and to recast the Einstein- Vlasov system into an autonomous finite- 
dimensional dynamical system over a compact state space. We note that the methods we use to 
prove our main result are closely connected with the general formalism developed in [5]. However, 
ensembles of massive collisionless particles do not directly fall into the class of matter models 
considered in [5] , which makes it necessary to generalize the approach. 

The paper is largely self-contained. In section [2] we discuss the collisionless matter model and 
give a derivation of the Einstein- Vlasov equations for the class of cosmological models we consider; 
the symmetry assumptions (LRS Bianchi type IX) are explained. In section [3] we reformulate the 
equations in terms of dimensionless variables; however, another reformulation of the equations is 
necessary to obtain a regular autonomous finite-dimensional dynamical system over a compact 
state space. The analysis of this dynamical system is performed in section [5] At the end of this 
section we state the main theorem: We prove that the a- and the w-limit set of generic orbits of 
the dynamical system that correspond to LRS Bianchi type IX solutions of the Einstein- Vlasov 
equations is a heteroclinic cycle. In the concluding remarks, section [5] we illuminate the main 
result from a physical perspective by putting it into a broader context. 



2 Derivation of the equations 

Consider an ensemble of massive particles that are in geodesic motion in a 'spacetime' (M, 4 g), 
which is a smooth four-dimensional manifold equipped with a metric tensor field of Lorentzian 
signature ( — h ++)• The assumption of geodesic motion reflects the condition of absence of in- 
teractions between the particles other than gravity; since the particles interact solely through the 
gravitational field they create collectively, the governing equations are Einstein's field equations of 
general relativity. This type of matter is commonly called 'collisionless matter', since, in particular, 
interactions by collisions are excluded. Examples of physical systems that are believed to be well 
approximated by the collisionless matter model in gravity are galaxies or galaxy clusters; in the 
former case, the particles are the stars of the galaxy, while in the cosmological setting, the particles 
are the galaxies of the cluster [1] |2] [15] . 

The ensemble of particles is represented by a 'particle distribution function' f p : R + x UM — > [0, oo), 
where UM C TM is the bundle of unit mass shells (four- velocity hyperboloids) , i.e., the subset 
of the tangent bundle TM given by the condition 4 g(it, u) = — 1 (where u is a future-directed 
four- velocity). Let x € M and u E UM be a four- velocity at x; let vo\um denote the induced 
volume element on UM; then 

/ p (m, x, u)vol UM dm 

represents the proper number density of those particles whose mass is in an interval of (infinitesimal) 
length dm around m and whose four-velocity is in an (infinitesimal) volume \o\um containing u. 
The proper (rest) mass density of the particles whose four- velocity is in a volume vo\jjm containing 
u is 

f(x, u) voIum = ( / mf p (m,x,u)dm) vo\ UM ; (1) 
this relation defines the 'mass distribution function' /. 
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Let (t,x l ) be a system of coordinates on M and {eo,et} be a frame, e.g., the coordinate frame 
{dt,d x i}; we assume that eo is (future-directed) timelike and spacelike, i = 1,2,3. Then the 
spatial components u l of the four-velocity (w.r.t. the frame) are coordinates on the hyperboloid UM 
and we can express the invariant measure on UM as voIum = y/\ dct 4 g| |uo| _1 du 1 du 2 du 3 , where 
Mo is determined from u l by the normalization relation 4 g(it,u) = g^ v u^u v — — 1. We adhere 
to the convention that spacetime indices are denoted by Greek letters, whose range is 0,1,2,3, 
while Latin indices are spatial indices and take the values 1, 2, 3. We use the Einstein summation 
convention. 

Regarding / (and / p ) as functions of the coordinates t, x 1 , u J (and to), i, j — 1,2,3, we find that 
the energy-momentum tensor of the ensemble is given by 

r"" = J vol UM J dmmf v u^u v = J vol UM f u^u" = J f u»u v yj\ det 4 g| \uo\~ 1 du 1 du 2 du 3 , (2) 

where u — w (w 1 ,m 2 ,m 3 ) by the condition 4 g(it, u) — g^ v u^u l/ = — 1. 

Remark. A common special case is the case where the ensemble of particles consists of one species 
of particles with equal mass m, which corresponds to / p (to, x, u) = 5(m— m) f (x, u). In this context 
it is customary to use the four-momentum v = mu as the variable of the distribution function; 
taking the relation between the volume elements on the unit mass shell UM and the mass shell 
PM = {v | 4 g(v, v) = — m 2 } into account, we see that f(x,v)vo\pM with f(x,v) — m~ 3 f(x, u) 
represents the proper number density. Then f(x,u) = m 4 f(a;, v) and ^ becomes 



f „r „\ ,, „ n—-, — ^—rdv 1 dv 2 dv 3 t. . „ „ r— — : — r dv 1 dv 2 dv 3 
J /(x,^)^V|det 4 g| = J f(x,v)v^v"y/\det* 



where vq is determined from v % by the mass shell relation 4 g(w,w) = g^uV^v" = — m 2 . 

Remark. The formalism to describe ensembles of massless particles is analogous. In the massless 
case, the domain of the distribution function is the bundle of future light cones, i.e., the set of 
future-directed null vectors. Accordingly, the energy- momentum tensor is given by ((2|), where Uq 
is determined from u 1 by the condition 4 g(u,it) = g^u^u" = 0. 

Both the particle distribution function / p and the mass distribution function / satisfy the Vlasov 
equation 

9t/+^/-V> V ^/ = > ( 3 ) 

which reflects the condition of geodesic motion of the particles; are the Christoffel symbols. 
Note in particular that the characteristic curves of the Vlasov equation, along which / is constant, 
coincide with the lift on UM of the spacetime geodesies. The gravitational interaction of the 
particles is modeled by the Einstein- Vlasov system, i.e., by coupling (|3]) to the Einstein equations 
of general relativity, 

Rfif — -^9nvR = T^ v . (4) 

In these equations, is the Ricci tensor of the metric 4 g and R = g^ v R^ v the Ricci scalar; for 
T^v we use the energy-momentum tensor ^ representing the Vlasov matter. We adopt units such 
that 8ttG = c = 1, where G is Newton's gravitational constant and c the speed of light. 

In this paper we consider spatially homogeneous spacetimes of Bianchi type IX that are locally 
rotationally symmetric (LRS), i.e., spacetimes of the form M = I x S 3 (where / is an interval of 
R) with metric 

4 g = -dt 2 + gil (t) lo 1 ® lj 1 + g 22 (t) (Cj 2 ® to 2 + w 3 <8> w 3 ) , (5) 

where {w 1 , w 2 , uj 3 } is a time-independent coframe on S 3 that satisfies dcu 1 = —uj 2 A ui 3 (and cyclic 
permutations) . As proved in [13] , the general solution of the Vlasov equation |3| on a background 
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spacetime with metric ([5]) can be expressed as 

/ = /oK,M 2 + M 2 ), (6) 

where ui = gijU J , i = 1, 2, 3, and fo : K x R + — » R + is an arbitrary sumciently smooth function. 

Let us compute from ([6]). Denoting by g the spatial Riemannian metric we have | det g| = det g 
and we find du 1 du 2 du 3 = (det g)^ 1 duidu 2 du 3 ; moreover, |i*o| 2 = l+g n u 2 +g 22 (u^+u 2 ). Therefore, 
the energy density p = T 00 is given by 

p = {detg)- 1 ' 2 J f (l+g 11 u 2 + g 22 (u 2 + u 2 )) 1/2 d Ul du 2 du 3 (7a) 
and the principal pressures p\ = T\ , p 2 — T\ , p 3 — T 3 3 are 

Pi = (dctg)- 1 ^ 2 f fog u u 2 (l + g^uj + g^iul + uDy 112 d Ul du 2 du 3 , (7b) 



where there is no summation over i. Since fo is the function ([6]), we find p 2 = p^. 

Remark. The energy density and the principal pressures ((TJ) depend on the arbitrary function 
Jo- This function can be interpreted as the 'initial data' for / at some time t — to, because 
f(to, u 1 , u 2 , u 3 ) = fo(gn(to)u 1 , (g22(io)) 2 ((w 2 ) 2 + (u 3 ) 2 )). Once the initial data fo is prescribed, p 
and p\, p2 are functions of the metric components, which can be interpreted as implicit relations 
between the principal pressures and the energy density, i.e., as an 'equation of state'. To compare, 
let us briefly recall the perfect fluid matter model. In the perfect fluid case, once an equation of 
state p = p(p) is prescribed, the energy density and the pressure are determined by the constraints. 
The evolution equations for the matter — the Euler equations — are equivalent to the conservation 
of the energy momentum tensor, V^T^^ = 0, and are thus contained in the Einstein equations (j4| 
(through the Bianchi identities). In contrast, the Vlasov equation ((3]) is independent; initial data 
for the first order equation (|3]) has to be prescribed in order to obtain a solution. 

For the pair ©-(J!]) to be a candidate for a solution of the Einstein- Vlasov system, the energy 
momentum tensor must be compatible with the LRS assumption, i.e., diagonal and T 2 2 = T 3 3 . 
This can be achieved by restricting to distribution functions (|6]) that are invariant under the 
transformation u\ — > —u\. These distribution functions are called reflection symmetric, see |16j . 
Besides reflection symmetry, for technical reasons, we also assume that fo has split support, which 
means that the support of fo does not intersect any of the axes. 

Remark. The (rest) mass current density of particles is given as 

= (det g)- 1/2 { fou^ \uq\~ 1 d Ul du 2 du 3 



and satisfies V ^N^ = 0, which expresses the conservation of (rest) mass. A straightforward 
consequence of the assumption of reflection symmetry is that N l — 0, i.e., the current density 
is orthogonal to the hypersurfaces t = const. The matter model can thus interpreted as being 
'non-tilted'. 

The Einstein equations (j4]) split into the Hamiltonian constraint equation 

9H 2 - ((fc 1 1 ) 2 +2(fc 2 2 ) 2 ) +R = 2p, (8a) 

and the evolution equations 

>h<f -'/••\.'/'. (8b) 
d t k\ = R\ - 3Hk\ - Pl + ~ (pi + 2 P2 ~ p) , (8c) 
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where there is no summation over i. The quantities k 1 1 and k 2 2 = k\ are the components of 
the second fundamental form of the hypersurfaces t — const, the Hubble scalar H is defined by 
H = -i trfc. 

1 (n 22 ) 2 1 (n 22 ) 2 

1 2 g 11 2 g 11 

are the non-zero components of the Ricci tensor, and R = R 1 1 + 2R 2 2 is the Ricci scalar. There 
is another constraint equation, the momentum constraint, that is satisfied identically by our as- 
sumptions. 

It is well known that the maximal interval of existence of solutions of the Einstein- Vlasov sys- 
tem (JT])-© is of the form (t_,i + ), where \t±\ < oo; we refer to [TTJUH]- After a time translation 
we may assume that the singularity in the past is at t- =0. The Einstein vacuum equations 
are recovered from (JTJ)— dHJ) by setting fo = 0. The asymptotic behavior of vacuum solutions is 
characterized by the existence of positive constants a±, b±, such that 

Sii (t) =a±(t-t±) 2 (l + o(l)) , S2 2 (t)=fl33(t)=&± + °(1) (*->*±). (9) 

This behavior is not exclusive to vacuum solutions; on the contrary, © is the typical behavior 
of solutions associated with a large variety of matter sources, (non-stiff) perfect fluids being the 
prime example [S]. Since the metric 

T: -dt 2 + a(t-t±) 2 (dx 1 ) 2 +b((dx 2 ) 2 + {dx 3 ) 2 ) (10a) 

with a,b — const is the so-called Taub solution (flat LRS Kasner solution), one refers to ([9]) as 
an approach to the Taub solution. (However, the spacetime associated with the metric (|10a|) is 
(0, oo) x T 3 instead of I x S 3 and thus of Bianchi type I; the Taub solution does not satisfy (|8d|).) 
The second class of vacuum LRS solutions of Bianchi type I (i.e., on (0,oo) x T 3 ) is the class of 
non-fiat LRS solutions 

Q: ~dt 2 + a{t-t±)- 2 ' 3 {dx 1 ) 2 + b{t-t ± f' 3 {{dx 2 ) 2 + {dx 3 ) 2 ) ; (10b) 

as opposed to (|10ap , the class Q does not play any particular role in the context of the asymptotic 
dynamics of vacuum solutions (or perfect fluid solutions) of (J?])-©. 

In this paper we are interested in the asymptotic behavior of solutions of the Einstein- Vlasov 
system ([7 |l — ([8]) as t — > t±. Our main result can be informally stated as follows: In the limit t — > t±. 
generic solutions of the system |7 ]) -([8 ]) oscillate between the Taub class and the class of non-flat 
LRS Kasner solutions. The asymptotic behavior of solutions of the Einstein- Vlasov system is thus 
qualitatively different from that of vacuum solutions. Therefore, "collisionless matter matters" , as 
opposed to (non-stiff) perfect fluid matter and several other matter sources which do not affect 
the structure of the singularity. 

Remark. The fact that the behavior toward the singularity of cosmological models with collisionless 
matter is qualitatively different from that of vacuum and perfect fluid models is known from LRS 
Bianchi type II models [TU [20] . The occurrence of oscillations in the asymptotic dynamics of LRS 
models that are induced by the anisotropy of the matter model has subsequently been studied in 
some detail, e.g., in 0[9]. 



3 Reformulated Einstein-Vlasov equations 

Let 
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Note that due to HI and ([6]), n = n(t) £ (0, oo) is proportional to the mass (or particle) density 
of the ensemble of particles (as measured w.r.t. the frame associated with (|SJ)). The variable s, 
on the other hand, satisfies s = s(t) 6 (0, |) and is a (non- linear) measure of the deviation of the 
metric from isotropy; s = -| (<*=> gn = 1722) corresponds to an isotropic metric. We further define 



1 



1 + 71 2 /3 



Since n 2 ^ 3 = (det g) -1 ^ 3 , n 2 ^ 3 corresponds to a length scale of the metric, and £ is a non-linear 
measure of such a scale; I = corresponds to a singularity (det g — 0); £ = 1 to a state of infinite 
volume (det 5 = +00). The equation satisfied by £ follows directly from (I8bp , 



= 2H£(l - t) 



(11) 



recall that H = — | trfc is the Hubble scalar; accordingly, H > means expansion, H < 
contraction. For s we find dts = — 2s(l — 2s)(fc 1 1 — fc 2 2 ), where /c^ — k 2 2 can be identified with 
(three times) the 2-2-component of the shear tensor. 

We are able to express the energy density (|7ap and the principal pressures (|7b[) in terms of £ and 
s; defining 

Pi p 1 pi+ 2p 2 1 



P 



w = - = - 

P 3 



= g (u>i + 2^2) 



(12) 



we obtain 



wi = - 2s) ■ 



/ /„ «? £(s 2 (l - 2s)) 1/3 + (1 - *)((! - 2s)u 2 + s(u 2 + «§)) 



-1/2 



rfu 



j /o 4 s2 (i - 2s )) 1/3 + (1 - - 2s K + -< u 2 + U D) 



1/2 



du 



w 2 = (l-£)s 



J f u 2 £(s 2 (l ~ 2s)) 1/3 + (1 - £)((! - 2s)u 2 + s(u 2 + u 2 )) 



-1/2 



du 



J f £(s 2 (l - 2s)) 1/3 + (1 - £)((! - 2s)u\ + s(u 2 + u 2 )) 



1/2 



(13a) 



(13b) 



du 



where du abbreviates du\du2du^,. The assumption of split support assures that u>i, W2 are smooth 
functions, since the denominator in (TIB")) is strictly positive for all values of £ and s. 

For our analysis it is necessary to recast the Einstein equations © into a different form. We use 
the dominant variable 



D 



H 2 



1 



3 J V 9 

to define normalized dimensionless variables according to 



(tr/c) 2 



1 



■9* 



(14a) 



Hr 



H 
D 



ft 1 ft 2 
3L> 



Mi 



1 5 



22 



f2 



3L> 2 



(14b) 



In addition we replace the cosmological time t by a rescaled time variable r via 



d_ 



Ddi 



(14c) 



Rewriting the Einstein equations (jHJ in the new variables, we obtain a decoupled ODE for D and 
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a system of ODEs that we call the reduced dynamical system: 



H' D = -(l-H D )(q-H D X + ), (15a) 

£' + = -(2 - q)H D ^+ - (1 - H 2 D )(1 -£+) + ! Mi 2 + « M*, «) - *)) , (15b) 

M{ = Mi - 4S+ + (1 - ££)£+) , (15c) 

f = 2H D £(1 - I) . (15d) 

In the system (fT5]l . g is the so-called deceleration parameter, q — 2£+ + |(1 + 3w)£l; in addition, 

f2 is determined from the variables £+ and Mi by the Hamiltonian constraint (|8a|) . and (fT4|) is 
used to express s as a function of Hd and Mi, i.e., 

O-l-fil-ijfl, .- (2 + 22^1)-'. (16) 



The system (fTS")) is a closed system that completely describes the dynamics of LRS Bianchi type IX 
Einstein- Vlasov models. Note that the r.h.s. of (fTS")) contains the functions (|13l) that are determined 
by an integration over fo = fo(ux, U2, U3). This means that, in particular, the r.h.s. of the dynamical 
system (|15l) depends on the initial data /o, which is an interesting feature of the problem. A more 
detailed derivation of (fT"5]) is given in [5j. 

Remark. The dynamical system (I15|) is invariant under the discrete symmetry 

t y t , H D -> -H D , £+ -> -E+ . 

Hence the qualitative behavior of solutions in the limit r — > +00 mirrors the behavior at r — >■ —00, 
and we may thus restrict ourself to study the latter. 

The state space for the dynamical system (|15|) is given by 
Six = Xix x (0, 1) , where X lx = {{H D> £+, M x ) | H D G (-1, 1) , M x > , Y? + + ^M\ < l} , 

see Fig. [TJ The set fix is relatively compact; the system (175]) is smooth on fix- However, (|T5"j) 
does not admit a regular extension to the entire boundary d£ix, which is because the variable s 
does not have a well-defined limit when Mi —> and Hp —> 1 simultaneously. This defect will be 
remedied by introducing the equivalent system (I19[) . 

Remark. Let us briefly comment on the equations describing the dynamics of cosmological models 
with different matter sources. The reduced dynamical system for perfect fluid matter is obtained 
from (|15|) by formally setting w\ — W2 = w = const in (|15b[) . which reflects the isotropy of the 
matter model (and the assumption of a linear equation of state) . The equation for i decouples from 
the remaining equations and the reduced dynamical system becomes the set of equations (|15ap ~ 
(|15c|) on the state space Aix- The reduced dynamical system in the case of an ensemble of massless 
particles is characterized by a decoupling of the equation for £ as well. This is because, in the 
massless case, the renormalized principal pressures are obtained from (|13[) by formally setting 
t = 0. The reduced dynamical system thus becomes the set of equations (|15a|) — (|15c|) . We therefore 
find that Ajx is the state space for both the perfect fluid case and the massless Vlasov case. We 
remark that this is not so for the lower Bianchi types. As shown in [TS], the state space for massless 
Vlasov particles has one dimension more than the state space for perfect fluids when the Bianchi 
type is I, II or III. In the Bianchi type IX case, the state spaces are identical, but another difference 
occurs: While the reduced dynamical system for perfect fluids admits a smooth extension to the 
boundary of ^ix, the Vlasov case is defective in this respect. Loosely speaking, the dynamics of 
Vlasov matter for massless particles does not live naturally in the state space Xjx- 
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Figure 1: The set X\x- 



4 Basic lemmas 

We use the system (fT5|) to prove two basic lemmas. 

Lemma 1. For every Bianchi type IX solution with Vlasov matter there exists tq € R such that 

• Hd{t) > Vt < To and Hd(t) is bounded away from zero as r — > — oo; 
. H d (t ) = 0; 

• Hd(t) < Vt > To and Hd{t) is bounded away from zero as r —> oo. 

Proof. The average pressure p of collisionless matter is non-negative; eqs. ([T2j) and (fT3|) imply 
w = p/p > 0. Therefore the general result by Lin and Wald [12] applies: Every Bianchi type IX 
model (with collisionless matter) possesses an initial singularity ('big bang'), expands initially, 
then reaches a time when the spatial volume is maximal, and finally recontracts to terminate in 
a singularity ('big crunch'). Hence there exists r such that Hd(t) > Vt < r , Hd(tq) = 0, 
and Ho{t) < Vr > To. It remains to prove that there exists a positive e such that Hrj(r) > e 
(Hd{t) < — e) for all sufficiently small (large) t. To do so assume the contrary, i.e., the existence 
of a solution whose a-limit set has a non-empty intersection with the plane Hd = 0. We first note 
that 

which is negative unless O = and £+ = 0. Second, 

|ff D =o,n=o,£ + =o = ) H D |Hij=o,n=o,s + =o = > H D |H D =o,n=o,s + =o = — 36 < . 

Let P be an a-limit point of the solution such that P lies on Hd = 0. Together with P, the entire 
orbit through P is contained in the a-limit set; hence there exists a sequence of times (t„)„ 6 n with 
t„ — > — oo (n — > oo) such that Hd{t,i — 5) > 0, i?rj(r n ) = 0, Hd{t u + S) < for a sufficiently small 
6] a contradiction. □ 

Remark. A posteriori, by Lemma U we find Hd(t) —> ±1 as r —¥ =poo. A note of caution: 
There exist anisotropic matter models other than collisionless matter such Lemma Q] is false; for 
such matter models there exist (typical) solutions that do not recollapse but expand forever (i.e., 
H d {t) > Vt); we refer to 0. 



Lemma 2. Let P be an a-limit point of a Bianchi type IX solution with Vlasov matter as repre- 
sented by an orbit of (|15p . Then 

V = o. 



Proof. The result is a simple consequence of eq. (|15dj) for ^ and Lemma [TJ □ 

Remark. Lemma [5] has a straightforward physical interpretation. Since (|15|) with £ = describes 
the dynamics of solutions of the Einstein- Vlasov equations with massless particles, see the remark at 
the end of section[3J we find that the asymptotic dynamics toward the past (and future) singularity 
in the massive Vlasov case are governed by the the dynamics of the massless Vlasov case. 



5 Analysis and main result 

By Lemma [TJ the subset Hd > of the state space Six is past-invariant under the flow of (|15P; 
this makes it possible to transform (|15|) to a different system of equations that is equivalent to (|15|) 
on the subset Hd > 0. 

Let 

Mf = 3rsintf, (17a) 

(l-#f,) = 2rcosi?, (17b) 

£+ = unchanged , (17c) 

where we assume that (Hd, Mi, £+) € X^ = X\x H {-Hd > 0}. The transformation of variables 
from (ifo, Mi, £_|_) to (r, £+), where r > and < i? < | is required, is a diffeomorphism on 
the domain (Hd,M\,H+) £ X^. We define the domain J^x °f t ne variables (r, £+) to be the 
preimage of that domain under the transformation (fl7|) . We obtain rcos ■& < \ from (IT7b)) ; the 
constraint 1 — — My = 17 > implies r sini? < 4(1 — £+)■ Therefore, J^ix can be written as 

y& = {r >0, € (0, S+ € (-1, 1) 1 r < min [— L^, ^ f +) ] }, (18) 



2 cos # sin $ 



see Fig. 2(a) 



In the new variables, the dynamical system (|15|) takes the form 

r' = 2r - H D £ + ) - 3S+ sin 2 4) , (19a) 

tf' = -3£+ sin(2tf) , (19b) 

= r sintf - 1 + (H D - S+) 2 + H D Y; + {q - H D Z+) + (iu 2 (4 «) - w x (l, s)) , (19c) 

i' = 2H D l(l-i) ■ (19d) 

where g = 2E^ + |(1 + 3iy)fi and 17 = 1 — — j r sin??. In addition, H D and s are regarded as 
functions of r and -d in (fHj|) . 

1 tan?? 



H D = Vl - 2rcos?? , s = - r . (19e) 

2 1 + tan?? v ' 

Note that these are smooth functions of r and i? on the state space J^x- In particular, 

ds 1 „ ,„ „ , 1 1 , 1 9s _ 1 , , . r „ 7T, 

— — = ; — - 2s 1 - 2s) = ; — r , hence - < — < - Vi? € 0, - 

dd sin(2i?) V ' 2 1 + sin(2i?) ' 4 ~ dti ~ 2 1 ' 2 J 
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(a) The space (b) A schematic depiction of 9Kt 

Figure 2: The space 3^rx an d its boundary. The various components of the boundary are defined 
in eqs. ([20]). 



and ds/dd = \ at = and d = § . 

The flow of the system (|19[) on the state space 3^ x (0, 1) is well defined in the past direction of 
time. This reflects the past invariance of the original domain x (0, 1). Accordingly, an orbit 
of (fT9|) represents the expanding phase of an LRS Bianchi type IX solution with Vlasov matter; 
conversely, the expanding phase of every LRS type IX model is represented by an orbit of (fT^j) . 

In contrast to the system (fT5j) . the dynamical system (fl9|) on the state space J^x x (0, 1) admits 
a regular extension to the boundaries. By the boundaries of (|18j) we mean 

{tfG (0,|),£+e (-1,1)}, (20a) 

{re(0,~),£+e(-l,l)}, (20b) 

{re (0,4(1 -£=)),£+€ (-1,1)}, (20c) 
1 4(1- £2 )i 4(1- £2 



r = 


: 


= 





= 


7T 






2 


n = 






2 cos sin $ 
and the closures of these invariant sets. The set 

1 4(1 -£ 2 + )- 



sini? 



,tfe(o,|),£+e(-i,i)}, (20d) 



2 cos -d sin <& 



2 cos 



^,*€(0,5),E + 6(-l,l)} 



is the preimage of the plane = under (fT7|) : it is not an invariant set; by Lemma [T] it is 
irrelevant for our considerations. This exhausts the list of boundaries of (fig)) . 

Lemma 3. Let P be an a-limit point of an orbit of (|19p (i.e., of a Bianchi type IX solution with 
Vlasov matter). Then 

r\ P =0 or tf|p = - . 



Proof. We consider the function 



^6 tan^> 
r 13 cos J v 
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which is smooth and positive on J^rx an d on the boundary subset $7 = 0, see (|20d|) . We find that 

Y' = -12Y(J? + + Q), Y'" ls+=QiQ=0 = -2iY(3 + H 2 D ), (21) 

hence Y is strictly monotonically decreasing along every orbit. This excludes the existence of 
a-limit points in J^x and on the boundary Q = 0. Suppose that there exists an orbit that 
possesses an a-limit point P with $| P = (and rip > 0). Then there exists a sequence of times 
(TVi)neN with r n — y — oo (n — ► oo) such that Y(t„) — > (n — > oo) along the orbit; a contradiction 
to (ED. □ 



Lemma [3] suggests to study the dynamical system that is induced by (|19[) on the boundaries r = 
and = | , see (|2TIaj) and (|2"0cl) . 

The boundary subset r = 0. The dynamical system p^|) induces the system 

i?' = -3S+ sin(2tf) , (22a) 
£V = -(i -£+)[£+- MM- ^2(0, s))] , (22b) 

on the boundary subset represented by r = and £ = 0, see (|20al) . From (Til?j) we see that 

, , , J/oU? [(l-2 S ) M 2 + s ( u 2 +u 2 ) l-l/2 rf 

wx(0,a = 1-25 1 LV — — , 23a 

ff [(l-2s)ul+s(u* + ui)} 1/2 du 

, s f/ uH(l-2 S ) U 2 + , s ( u 2 +u 2 ) l-l/2 d 

w 2 0, s = a J - 2 11 — — %- 9 , 23b 

Jf [{l-2s)ul + s(ul + u 2 3 )] 1/2 du 



which implies that w — + 2w 2 ) = \- Recall that in the context of (j2"2")l . s = s{d), see 

It is straightforward to see that iui(0, s) — u>2(0,s) = 1 — 3^2(0, s) is a strictly monotonic func- 
tion [T5]. Since ^(O, 0) = and 1^2(0,5) = |, there is a unique value s G (0, 5) such that 
^2(0, So) = \- Let ^0 denote the (unique) value of & associated with s by (|19e[) . We conclude 
that there exists a unique fixed point F in the interior of the set r = 0, £ = 0: 

F The fixed point F is given by r = 0, i9 = # , = 0, £ = 0. 

The remaining fixed points of the system (|2"2"|) are located on the boundary of the subset r = 0; 
these are given by # = 0, $ = 5, and f2 = (with E + = ±1). 

Tj The Taub point Tj is given by r = 0, i9 = ^, = — 1, £ = 0. 

Qj The non-flat LRS point Qj is given by r = 0, $ — -|, = 1, £ = 0. 

Rjt The point R s is given by r = 0, = f , £+ = |, £ = 0. 

T b The Tawfe pomt T b is given by r = 0, ■& = 0, S + = -1, £ = 0. 

Qb The non-flat LRS point Q|, is given by r = 0, = 0, £+ = 1, £ = 0. 

The physical interpretation of the fixed points is straightforward. At the fixed point F, since 
w\ = W2 — w 3 = w — the principal pressures are identical and the matter is thus isotropic. 
Since E + = 0, we have k\ = k 2 2 = k 3 3 , i.e., the geometry is isotropic as well. Accordingly, the 
fixed point F represents the Friedmann-Lemaitrc- Robertson- Walker model, 

3n = 322 = 533 = at (a > 0) , (24) 



12 



O 



Qb 



Figure 3: The flow induced on the two-dimensional boundary subset determined by r — and 
I = 0. The color-coding of the fixed points has the following meaning: A fixed point depicted in 
black (resp. white) is an attractor (resp. repeller) of orbits in a transversal direction (within the 
set I = 0). 

i.e., a spatially flat and isotropic solution of the Einstein- Vlasov system. 

The fixed points T\, and Tj, as well as the orbit connecting these two fixed points, cf. Fig. [3J 
represent solutions of the Taub class, i.e., (jlOal) . During the time that an orbit of (p~5|) is close to 
T|,, Tjj, or Tb — > Tjj, it represents a metric (|5|) whose components are well approximated by (|10a|l . 
(To see this we simply use the variable transformations (fT^j) and ((T7)) .) Likewise, while an orbit 
of (fl~5)) is close to Q|,, Qj, or Q|, Qj, it represents a metric ([5]) whose components are well 
approximated by (I10b[) . 

The solution associated to the fixed point Rj is of type I as well. (This fixed point has already 
been found in [T!5], where it was denoted by P3.) The metric is 

,9ii = a , 522 = 533 = bt 4/3 (a, b > 0) ; 

we refer to [5]. Note that this a non- vacuum solution and that the matter is anisotropic; hence Rj 
does not correspond to a perfect fluid solution. 

Lemma 4. The system (|22[) on the closure of the subset r — (and t = 0) gives rise to the flow 
depicted in Fig. [3J 

Proof. Consider the function 

Z = (1 - S^r 1 ((1 - 2s) s 2 )) - 1/6 1 /o [(1 - 2s)u{ + s{ul + ul)] 1/2 dv , 

which satisfies 

Z' = -1ZY? + , ^'|'e+=o = -36Z(w 2 (0,s) - i) 2 . 

Since Z is strictly monotonically decreasing except at F, we conclude that F is the w-limit set 
of every interior orbit, while the a- limit set is contained on the boundary. The straightforward 
analysis of the flow on the one-dimensional boundaries leads to the claim of the lemma. □ 

The boundary subset i9 = ^. The dynamical system (|19j) induces the system 

r' = 2r(Y? + - 4E+ + 1 - I) , (25a) 

K = 5(2 - E+) - (1 - X\ - j) (E + - ~) , (25b) 
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Figure 4: The flow induced by ([T9|) on the two-dimensional boundary subset determined by $ = § 
and £ = 0. 

on the boundary subset represented by $ = \ and ^ = 0, see (|20c[) . Here we have used that s = \ 
when § = see (|19el) . which implies that w»i(0, s) = and W2(0, s) = |, see (I2U1) . 

The one-dimensional boundary of the system (f2l)j) consists of a part where O = and a part r = 0; 
on the latter, the systems (|22|) and (|25l) intersect. 



Lemma 5. The system (|25|) on </ie closure of the subset •d = \ (and £ = 0) gives rise to the flow 
depicted in Fig. [7} 

Proof. The absence of fixed points, periodic orbits and heteroclinic cycles in the interior of this 
set implies, by the Poincare-Bendixson theorem |14) . that the a- and w-limit points of orbits must 
be located on the boundary. A local dynamical systems analysis then yields the claim of the 
lemma. □ 

This completes our analysis of the boundary subsets. Let us return to the study of (|T9|) . 

Lemma 6. There exists a one-parameter family of orbits of (fTi?|) whose a-limit set is the fixed 
point F. 

Proof. A simple calculation shows that the unstable manifold of F is two-dimensional. (See 
Lemma |4] for the stable manifold.) □ 

Lemma 7. Assume that a point P of the (interior of the) boundary subset r = 0, I = 0, see (|20a|) . 
is an a-limit point of an orbit of (|19j) . Then P = F. 



Proof. Let the orbit under consideration be denoted by 7. Assume that P 7^ F. Together with P, 
the entire orbit through P and its cj-limit point F must be in the a-limit set of 7; see Lemma HJ 
Since F is a saddle point, see Lemmas [4] and [6j there is a point P on the unstable manifold of F that 
is contained in ad). (Lemma [2] enforces P to be located on I = 0.) Let (r n ) n6 N with r„ — > —00 
(n — > 00) denote a sequence of times such that 7(t„) — > P (n — > 00); by assumption, none of 
the points 7(t„) are contained on the unstable manifold of F. The orbit through P represents a 
solution of the Einstein- Vlasov equations that satisfies Lemma [TJ By continuous dependence on 
initial data we can thus construct a sequence of times f n with f n — > — 00 (n — > 00) such that 
Hd{tu) — ► (n — > 00); this is in contradiction to Lemma [TJ □ 

Lemma 8. There does not exist any point in the (interior of the) boundary subset $ = 5, £ = 0, 
see (I20c[) . that is an a-limit point of an orbit of ([T! 
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Proof. Assume that there exists an orbit 7 such that P € 0(7), where P lies on the (interior of 
the) boundary subset # = §, see (|20c|) . Together with P, the entire orbit through P and its w-limit 
point R(j must be in the a-limit set of 7; see Lemma [5] The fixed point Rj is a saddle point, 
see Lemmas 2] and hence there is a point P on the unstable manifold of R« that is contained in 
0(7). The unstable manifold of Rj is two-dimensional; a particular orbit contained in it is the orbit 
Rj — >• F depicted in Fig. [3] the orthogonal unstable direction is the ^-direction. Lemma [2] implies 
that £^p = 0, hence P is located on the orbit Rj — » F of Fig. [3] However, this is a contradiction to 
Lemma [71 □ 

Lemma 9. There does not exist any point on the (interior of the) line Tjj — > Rj Qj, see figs. [3] 
andJ^J £/ia£ is an a-limit point of an orbit of (|19[) . 

Proof. The proof is analogous to the proof of Lemma [H □ 

Summarizing the statements of the previous lemmas we obtain the main theorem of this work. We 
give a formal statement of the theorem using the developed framework; the less formal interpreta- 
tion of the theorem is presented in the concluding remarks. 

Theorem 1. Consider an orbit of (|19[) representing (the expanding phase of) a Bianchi type IX 
solution with Vlasov matter. Then the a-limit set of this orbit is either the fixed point F, which 
is the non-generic case, or it is the heteroclinic cycle of orbits connecting the four fixed points T^, 
Tj, Qtj, Qi,, see Fig. 0' this is the generic case. 

Proof. By Lemma [6] there exists a non-generic family of orbits whose a-limit set is F. Consider an 
orbit that is not a member of this family. Lemma [2] implies that the a-limit set must be located 
on the boundary I = of the state space of (fT9|) . Lemma [3] enforces a-limit points to lie on the 
closures of the boundary subsets r — and d — ^. However, Lemmas [7J and [5] imply that points in 
the interior of these subsets are excluded. Finally, by Lemma |H] we find that the set that remains 
as the possible location of a-limit points is a one-dimensional set, the heteroclinic cycle connecting 
the fixed points T^, Tj, Qjj, Qt,. □ 



6 Concluding remarks 

Let us give an interpretation of the main result, Theorem[TJ in informal terms: Every LRS Bianchi 
type IX solution with collisionless matter (i.e., solution of the Einstein- Vlasov system) possesses 
an initial singularity (which is chosen to be at t = 0). The behavior as t — > of generic solutions, 
where by 'generic solutions' we mean a family of solutions that corresponds to a set of initial data 
of full measure, is characterized by oscillations between the Taub solution (|10a[) and the non-flat 
LRS solution (|10b|) : There exists an infinite sequence of (increasingly small) time intervals such 
that the components of the metric ([5]) are well approximated by (|10a|) . and a sequence of time 
intervals where (|10bl) yields a good approximation; the accuracy of the approximation increases as 
(->0. During the transitions between (|10a[) and (flObJ the metric takes an entirely different form; 
we merely note that each transition from the non-flat LRS solution to the Taub solution (which 
correspond to an orbit of (fl9| closely following the orbit Qt, — > T^ in Fig. [5]) is characterized by the 
fact that the variable il attains values close to 1. This means that the influence of the matter on 
the (asymptotic) dynamics cannot be neglected. The oscillatory behavior toward the singularity 
of generic LRS Bianchi type IX cosmological models with collisionless matter is thus qualitatively 
different from that of perfect fluid models with the same symmetry. In fact, in the LRS case, perfect 
fluid cosmologies of Bianchi type IX are asymptotic to the Taub solution (|10a|l . which is also the 
behavior of LRS Bianchi type IX vacuum cosmological models. Therefore, as the 'structure' of 
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Figure 5: The projection of a generic orbit (blue line) of the dynamical system fp~9|) to the set 
I — 0. Dashed red lines are orbits that lie on the boundary of the state space and connect to form 
a heteroclinic cycle. There exist non-generic orbits that are asymptotic to the fixed point F (not 
depicted). 

the singularity is concerned, the assumption of perfect fluid matter does not change the vacuum 
behavior ('matter does not matter'), while 'collisionless matter matters'; Vlasov matter has an 
important effect on the structure of the singularity. 

The occurrence of oscillatory dynamics toward the singularity is intimately connected with the 
(in)stability of the Kasner fixed points (i.e., the Taub points and the non-flat LRS points in the 
context of the present work). While stable in the context of the simplest cosmological models, 
which are characterized by a small number of (true) degrees of freedom and by the isotropy of the 
matter model, the Kasner points run the risk of loosing their stability when degrees of freedom are 
added to the problem or when the matter model becomes anisotropic. Loss of stability is then a 
possible cause for the existence of heteroclinic structures which in turn induce oscillatory behavior 
of cosmological models. The prime example is the transition from (non-LRS) Bianchi type VIo 
and VIIo models to (non-LRS) type VIII and IX models in the vacuum case. The dimension 
of the state space increases, the Kasner points lose their stability properties, and the dynamics 
of (vacuum) cosmological models become the oscillatory "Mixmaster" dynamics. Similarly in 
spirit, in the context of LRS Bianchi models with Vlasov matter, the dimension of the state space 
increases when one goes over from LRS type I to LRS type II models. The loss of stability induces 
oscillatory behavior of LRS Bianchi type II solutions with Vlasov matter, see [191 120) . However, 
already a change of matter model alone (where the number of degrees of freedom of the problem 
are unaffected) can change the stability properties of the Kasner points: For instance, while the 
asymptotic dynamics of LRS Bianchi type I solution with Vlasov matter are "monotone" , one 
observes oscillatory behavior for LRS Bianchi type I solutions with elastic matter [9] . The present 
work provides another example: The state spaces and reduced dynamical systems that represent 
the dynamics of LRS type IX perfect fluid models and LRS type IX models with massless Vlasov 
matter are of the same dimension; the same is true for the state space of LRS type IX perfect 
fluid models with non-linear equations of state and the state space of type IX models with massive 
Vlasov matter considered in this paper. Despite this fact, the qualitative dynamics are radically 
different, which is ultimately because the Kasner points exhibit different stability properties. In 
brief: "Matter matters" . 

Note that by the symmetry of the problem, the behavior of generic solutions toward the final 
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singularity (see Lemma Q} mirrors the behavior toward the initial singularity. Furthermore, we 
note there exists a non-generic family of LRS Bianchi type IX solutions with collisionless matter 
that isotropize toward the initial singularity, see Lemma [6] These solutions are asymptotic to the 
isotropic Friedmann-Lemaitre-Robertson- Walker model (|24)l . i.e., gij(t) ~ tSij as t — > 0. Although 
these non-generic solutions do not influence the asymptotic behavior of generic ones, they may 
affect the intermediate behavior of generic solutions and thus play a relevant role in physics. We 
refer to |21j for a discussion on this interesting topic. 

In this paper we have shown that collisionless matter has an important effect on the structure of 
the singularity already in the context of solutions of the Einstein equations with high symmetries 
(LRS Bianchi type IX). It is to be expected that this type of effect becomes even more pronounced 
when the degree of symmetry is reduced. In particular, an exciting question is in which manner 
collisionless matter will influence the so-called Mixmaster behavior of generic (non-LRS) vacuum 
Bianchi type IX models. Naively, one might expect that the two types of oscillations, the Mixmaster 
oscillations (which are induced, in a manner of speaking, by gravity itself) and the oscillations 
caused by the anisotropy of the matter model 'intertwine' to yield a intricate oscillatory structure. 
Unfortunately, it is difficult to give a well formulated conjecture, the reason being that it is not 
known at present how to extend the dynamical systems method to the full (non-LRS) Bianchi 
type IX case. 
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